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Elegance of disordered granular packings: a validation of Edwards' hypothesis 
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We have found a way to analyze Edwards' density of states for static granular packings in the 
special case of round, rigid, frictionless grains assuming constant coordination number. It obtains 
the most entropic density of single grain states, which predicts several observables including the 
distribution of contact forces. We compare these results against empirical data obtained in dynamic 
simulations of granular packings. The agreement is quite good, helping validate the use of statistical 
mechanics methods in granular physics. The differences between theory and empirics are mainly 
related to the coordination number, and when the empirical data are sorted by that number we 
obtain several insights that suggest an underlying elegance in the density of states. 

PACS numbers: 45.70.Cc, 05.20.Gg, 05.10.Ln, 05.65.+b 



The intriguing behaviors of sand and other granular 
materials are not well understood from a fundamental 
point of view [ij and there is no theory with a pedi- 
gree equivalent to the Navier-Stokes or Maxwell's equa- 
tions to explain how their state evolves over time, even 
in the most commonly occurring scenarios. Consider- 
ing the ubiquity of granular materials in nature, this is 
quite surprising. Making a new effort to explain their 
physics, Edwards and Oakeshott proposed that the meth- 
ods of statistical mechanics may be successfully applied 
2] . They hypothesized a priori a flat measure in the sta- 
tistical ensemble — that every metastable arrangement of 
grains (a blocked state) is equally probable under com- 
mon conditions — and that the analysis of this ensemble 
should predict some of the important behaviors. 
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Because dynamics of granular materials are nonlinear, 
lossy, and quite different than the dynamics of atoms, 
it is an important question whether they are ergodic or 
whether they bias the measure such that Edwards' hy- 
pothesis would not be correct. Seeking to answer this, a 
number of empirical tests have been performed by com- 
puter simulation. In these, Edwards' hypothesis has suc- 
cessfully predicted packing behaviors for several idealized 
models 3] and the diffusion-mobility behavior of individ- 
ual grains when a simulated packing is slowly sheared Q . 
In both cases it appears that the dynamics cause the ge- 
ometry of the model to explore some region of the phase 
space with sufficient ergodicity to justify the flat mea- 
sure. Also, experiments vibrating a powder have shown 
that it achieves a steady-state volume, repeatable but 
dependant on the frequency and amplitude of the vibra- 
tion Edwards has used that as the starting point 
to develop a Boltzmann equation 0, which assumes the 
individual grains occupy volumes of space that are sta- 
tistically uncorrelated to that of their neighbors. Surely 



FIG. 1: (Somilogarithmic) Distribution of granular contact 
force magnitudes Pj (/) (main graph) and Cartesian compo- 
nents Px(fx) (inset). The theory and empirical discrete ele- 
ment model (DEM) are strikingly in agreement. This implies 
that Edwards' hypothesis is sufficient to capture important 
organizational features of quasi-static granular physics. 



for friction-dominated packings (such as powders) this is 
reasonable, and so Edwards' transport equation proves 
ergodicity in the Boltzmannian sense. 

In this Letter we present a different kind of test for Ed- 
wards' hypothesis. Rather than examining the geometric 
features of the packing, we demonstrate that Edwards' 
flat measure correctly predicts the distributions for sin- 
gle grain load states and for contact forces. This pre- 
diction is centrally important to a statistical mechanics 
theory because the distribution reflects how the ensemble 
is organized and demonstrates whether or not the correct 
physics have been incorporated. In particular, Edwards' 
hypothesis should predict at least three features in the 
contact force distribution Pf(f) as illustrated in Fig. ^ 
the wide tail Q related to the heterogeneity of stresses in 
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a packing (force chains) [8(; the small peak near the av- 
erage value of force under isotropic conditions [9| related 
to static equilibrium of the grains (jamming) [l0(] ; and 
the non-zero probability density at zero force, P/(0) > 0, 
related to the tipping of grains (fragility) The com- 
bination of these three features is unique to the granular 
distribution, not being found in the typical densities of 
states for thermal systems. If Edwards' hypothesis fails 
to produce this form then it is doubtful that it could be- 
come the basis for a theory of quasi-static rheology, since 
the tipping or sliding of individual grains depends upon 
the state of their contact forces. 

Following Edwards and coworkers ^3] , we focus on the 
case of amorphous packings of cohesionless, rigid grains 
all having the same coordination number Z that makes 
the packing isostatic 0] ■ 0m case is further idealized by 
using two dimensional round grains with monodisperse 
grain diameters, omitting gravity and working in the 
thermodynamic limit (infinitely large packings) so that 
the boundary layer may be neglected. We focus on the 
frictionless case so that Z = 4, and we limit this Letter to 
isotropic stress and fabric although our methodology can 
solve for anisotropic cases, too. The idealizations may be 
taken out in future refinements of the theory, but this is 
a good starting point because packings of cohesionless, 
round grains that are perfectly rigi d Il 4| and /or friction- 
less and/or monodisperse jil ll(:il| have been the 
focus of many empirical studies and are known to have 
force distributions with the same features as the less ide- 
alized packings. Hence they must be subject to the same 
basic organization in the physics. 

The goal of the analysis is to combine Edwards' micro- 
canonical DOS and contact force probability func- 
tional and then derive the density of single grain 
states p g (w x , w y , Q\, . . . , #4). The first two arguments of 
this density, the Cartesian loads, are 

1 4 1 4 

Wx = 2 zZ^ cos9 f^ w y = 2 z2^ sin9 ^- W 

7=1 7=1 

where / 7 and 9 1 are the four contact force magnitudes 
and contact angles on a grain. In the special case we have 
selected, solving for p g provides everything that can be 
known about the individual grains in the packing. For ex- 
ample, it contains the joint distribution of contact forces 
and angles, 

/>OC />27T -1 4 

P fe (f,9) = d 2 w d 4 6 p e x-Y / S(9-6 7 ) 
Jo Jo 4 7=1 

x s[f - f 1 (w x ,w y ,e u ...,e 4 )}, (2) 

and the fabric distribution P4g(9\, . . . , 84) (lif . 

The analytical method 0] is to count states in Ed- 
wards' ensemble and maximize entropy applying the 
same well-known procedure that has been used to derive 




FIG. 2: Portion of a Discrete Element Model (DEM) showing 
the disordered packing fabric and propagation of force chains. 
Line width is proportional to force magnitude. Although dis- 
ordered, a simple pattern can be discerned in the density of 
states as discussed in the text. 

the Bose or Fermi distributions [2(J. The result is 

PgK, w v , 9 1 ,...,e 4 ) = G{6i, ...,9 4 ) e ~ x * w *- x y w y 
4 

^n[ p /»(^M i/29 (/))' ( 3 ) 

7=1 

where O is the Heaviside (unit step) function, and X y 
are the Lagrange multipliers that scale mechanical load- 
ing, and G derives from the array of Lagrange multipliers 
used to conserve P40. Eqs. © and © form a recursion 
in Pfg and p s , the "transport" equation, which may be 
solved numerically using P40 and the mechanical loading 
as inputs. 

For the present the transport equation has been solved 
in the isotropic case with a simplifying approximation: 

Pg(w x ,w y ,9f}) w p w (w x ,w y )pe(9f)) 

xes(w x ,Wy,6f)), (4) 

where 0s is a function that evaluates either to unity 
or zero if the grain is stable or unstable, respectively 
[2l|. This modified separability assumes no correlation 
between the loads and fabric apart from the truncating 
effect of 65. The physical idea is that correlation does 
arise predominantly because nature disallows unstable 
grains. Empirical results have shown this to be correct 
|l9| . In the remainder of this Letter, "the theory" refers 
to the resulting numerical solution. 

To compare with the theory, we have performed dis- 
crete element modeling (DEM) of 17,000 two dimen- 
sional, round, frictionless grains. A portion of our DEM 
packing is shown in Fig.|3to contrast the spatial disorder 
of its force network with the simplicity of the statistics, 
discussed below. The grain diameters are uniformly dis- 
tributed from 1.0 to 1.5 to reduce crystallization. The 
grains were deposited isotropically into a square test cell 
with hard walls and without gravity, and their diameters 
increased by rescaling, producing the desired isotropic 
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FIG. 3: Distribution of s = (w x — w y )/t related to the 
shear stress borne by single grain states. Solid curve — DEM. 
Dashed curve — theoretical prediction, the numerical solution 
of the transport equation. 

stress state. The grains were allowed to move dynami- 
cally until they located and settled into a blocked state. 
They have a linear spring contact law, but staying near 
the jamming transition avoids excessive deformation of 
the contacts and approximates the grain rigidity of the 
analysis. Data from grains in the boundary layer (cho- 
sen to be 4 grain diameters along each wall) were dis- 
carded to reduce the boundary effects, which we found 
will otherwise significantly skew the statistics. We also 
note that the theory assumes Z — 4 for every grain, but 
the DEM produces significant populations for Z = 3, 4, 
and 5. Therefore, to evaluate Edwards' hypothesis we 
will discuss the differences between these populations. 

Fig. ^ shows the DEM data compared to the theory 
for Pf(f) and for the distribution of Cartesian compo- 
nents of force, P x (fx)- They are in remarkable agree- 
ment, demonstrating all the correct features and thereby 
indicating that the ensemble naturally incorporates the 
correct contact force physics. 

To investigate the DOS more fully we note that w x 
and w y are not statistically independent and therefore 
we would need to plot their statistics as a joint distribu- 
tion. However, the change of variables to t ~ w x + w y 
and s = {w x — w y )/t achieves (approximate) statisti- 
cal independence so that they can be separated more 
meaningfully. The parameter t is analogous to hydro- 
static pressure but at the grain scale whereas s is a ra- 
tio that indicates the degree of shear stress at the grain 
scale. The distribution of the latter, P s (s), in Fig. |3 
demonstrates remarkable agreement between the DEM 
data and the theory. We can fit them to a functional 
form, P s (s) = cos(7rs/2)exp(-s 2 /2cr 2 ) with a = 1/4. To 
explore the dependence on Z we segregated the DEM 
data into Z = 3, 4, and 5 populations and plotted P s (s) 
for each in Fig. 21 Remarkably, a good fit to each popu- 
lation is made simply by writing a — 1/Z. This identifies 
a previously unknown pattern in the form of the DOS. 

This result has an interesting relationship with recent 
work on the statistics of cooperative bridges. These 
bridges naturally occur within the bulk of packings and 
have been a focus of much interest due to the way that 
they direct the propagation of stress. In their recent 
work, Mehta et al. found that the lengths of these 
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FIG. 4: The same distribution s from the DEM as in Fig. 
except segregated by grain coordination number Z into three 
graphs (shifted vertically for clarity). An empirical fit was 
suggested by the theory (dashed curves) which fits all three 
Z populations when the standard deviation a = 1/Z is the 
only parameter. 

bridges have an exponential distribution like P/(/), and 
so bridges were proposed to be a geometrical analog of 
the force chains themselves (2^. At the same time, they 
found the spatial orientations of the bridges to have a 
Gaussian distribution, which is similar to P s {s). It has 
been pointed out to us that this bridge orientation is in 
fact closely related to s because the angle with respect 
to the gravity vector determines the shear stress borne 
by a bridge. Thus, it appears that both the exponential 
and the Gaussian statistics found within the single grain 
DOS may be connected, through the mesoscopic feature 
of bridges, to important macroscopic behaviors. 

The other statistically independent variable, i, was also 
analyzed in the theory and found to have the distribution 

P t (t) = t^e-P* (5) 

where t has been normalized and where [3 — 5. This 
is an extremely interesting form when several facts are 
considered. First, it is well known in probability theory 
that, when several independent random variables U hav- 
ing distributions Pi are added together, T — ^ ti, then 
the distribution of the sum is 

P T =Pl®P 2 ®---®Pn (6) 

where ® is the convolution operator. Second, it is well 
known that Eq.JSJis the convolution of pure exponentials, 

£/3-i e -/3t = e -* ( g )e -*(g,...(g) e - t ) t>0 (7) 

where there are (3 exponentials being convolved, to be 
precise. Third, a pure exponential is of course the canoni- 
cal (Gibbs) distribution. Together, these facts tell us that 
the hydrostatic load on a grain in a disordered packing is 
distributed as if it were composed of several statistically 
independent, canonical contributions. This is quite sur- 
prising because the contact forces themselves are neither 
independent nor canonical. 

To check this, we segregated the DEM data by Z and 
obtained Pt{t) for each population as shown in Fig. [S] 
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FIG. 5: Distribution of t (hydrostatic loading of the grains) 
from the DEM, segregated into three graphs by their coordi- 
nation number Z (shifted vertically for clarity). The empiri- 
cal fits are Pt(t) = t l3 ~ 1 e~ ,3t as predicted by the theory, using 
P = 2Z - 4. 



This confirmed the pattern: all three populations are 
fit perfectly by Eq. [SJ using (3 = 2Z — 4 as the only 
parameter. While the origin of the value of (3 is yet to be 
explained, it is clear that the essential physics have been 
correctly incorporated into this theory because the forms 
of the distributions are all correctly predicted. 

Furthermore, an important feature of (3 can be seen: 
averaging over all the grains in the packing, (f3) = (Z) if 
and only if (Z) = 4. This also happens to be the con- 
dition for mechanical isostaticity and recent studies have 
demonstrated that it really is satisfied for the present 
case This means that, if the independent canonical 
variables suggested by Eq.UJcan be identified, we will find 
that the number of them is exactly equal to the number of 
contact forces in the packing. This is surprising because 
in general (3 ^ Z and therefore the new variables cannot 
be localized to the individual contacts. This forms an in- 
teresting analogy to the molecular vibrations in a solid, 
which are resolvable into non- localized phonon statistics, 
or to the eigenmodes of a mass-spring network. 

In summary, the theory predicts the correct forms for 
Pf(f), Px(fx), Ps(s) and P t (t). By segregating the grains 
of a dynamic simulation by their coordination number Z , 
we discover that all the populations fit the theory's pre- 
dicted DOS (represented by s and t) with only a simple 
parameter change based on Z. This identifies a previ- 
ously unrecognized but elegant pattern. 

We conclude that all of the features of Pf(f) (as pro- 
duced by dynamic simulations and experiments) are nat- 
urally predicted by Edwards' hypothesis, alone; none of 
these features are the result of dynamically-induced de- 
partures from a flat measure. Therefore, Edwards' hy- 
pothesis, without recourse to the individual grains' dy- 
namics, produces an ensemble that contains the force 
chains, the fragility, and all the other important granular 
phenomena that have been correlated to those features 
of Pf(f). The results of this Letter therefore provide one 
more indication that Edwards' hypothesis may be the 



correct starting point for a complete statistical mechan- 
ics theory of granular packings. 
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Bhattacharya of the University of Central Florida Physics 
Department and with Robert Youngquist of NASA's 
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